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Abstract 

J 

The  finite  element  method  was  used  to  solve  the  flow  field  problem 
around  a  thin,  lifting,  flat  plate,  airfoil.  The  governing  equation 
solved  is  the  laplace  equation,  which  is  valid  for  inviscid,  irrotalional, 
incompressible  flow. 

The  finite  element  equations  were  derived  through  the  method  of 
weighted  residuals  with  weighting  functions  selected  by  the  Galerkin 
method.  For  the  purposes  of  analysis,  the  infinite  flow  field  was 
replaced  by  a  finite  domain.  Neumann  type  boundary  conuitions  were 
imposed  on  the  airfoil  surface.  Dirichlet  boundary  conditions  were 
specified  as  required  by  the  problem  formulation  for  uniqueness. 

Three  types  of  solution  methods  were  used,  for  various  treatments 
of  the  jump  discontinuity  required  in  the  lifting  problem.  The  first 
method  was  a  superposition  technique,  which  treated  the  potential  along 
the  upper  and  lower  nodes  of  the  branch  cut  as  constant.  The  circulation 
was  determined  by  applying  the  Kutta  condition  during  the  combination 
of  the  subproblems.  The  second  method  was  an  iterative  technique  where 
the  circulation  was  varied  until  the  Kutta  condition  was  satisfied.  This 
method  also  specified  constant  potential  along  the  upper  and  lower  branch 
cut  nodes.  The  third  method  was  also  an  iterative  technique  on  circula¬ 
tion;  however,  only  the  ratio  of  potentials  across  the  branch  cut  nodes 

were  kept  constant.  - 

\ 

Three  types  of  elements  were  investigated  to  approximate  the  solution 
for  the  velocity  potential  function.  The  first  was  a  bilinear  rectangular 


v  i 


element.  The  second  was  a  mixed  element  with  three  linear  sides  and  one 
quadratic  side  used  only  on  the  airfoil  surface.  All  other  elements  used 
with  it  were  bilinear .  The  third  element  used  was  a  biquadratic,  Lagrange, 
element.  The  convergence  characteristics  of  each  element  were  studied 
as  a  function  of  the  discretisation.  The  pressure  distributions  are 
compared  with  those  of  classical  thin-airfoil  theory. 


FINITE  ELEMENT  ANALYSIS  OF  SUBSONIC 
FLOW  OVER  A  LIFTING  THIN  AIRFOIL 

I  Introduction 

In  the  design  of  aircraft  lifting  surfaces,  it  is  important  to 
select  a  design  that  has  a  high  probability  of  meeting  desired  charac¬ 
teristics.  Because  of  the  large  costs  associated  with  the  construction 
and  test  of  models,  computational  techniques  for  predicting  the  perform¬ 
ance  of  these  designs  have  evolved  with  aircraft  technology.  Over  the 
past  25  years,  the  most  extensively  used  method  in  the  aircraft  indus¬ 
try  has  been  the  finite  difference  method.  This  method;  however, 
becomes  difficult  to  apply  when  complex  geometries,  multiply  connected 
domains,  or  complex  boundary  conditions  are  involved.  Finite  element 
methods  can  overcome  some  of  these  difficulties,  through  easier  treat¬ 
ment  of  complex  geometries  and  a  more  consistant  method  of  using  higher- 
order  approximations  (Ref  l).  For  these  reasons  the  finite  element 
has  been  increasingly  utilized  in  solution  of  fluids  problems.  The 
treatment  of  the  flowfield  as  a  two-dimensional  potential  flow  problem, 
although  the  flow  is  in  fact  more  complex,  is  a  reasonable  simplification 
since  the  selection  of  a  suitable  section  shape  is  an  important  part  of 
the  design  process.  The  application  of  this  method  to  the  thin  lifting 
airfoil  is  a  good  first  step  in  the  solution  of  the  general  airfoil 
problem.  A  further  simplification  can  be  made  by  assuming  incompressible 
flow,  thus  reducing  the  problem  to  solution  of  the  Laplace  equation. 
History  and  Previous  Work 

The  finite  element  method  was  developed  by  aircraft  structural 
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engineers  in  the  1950's,  to  analyze  large  structural  problems  in  air¬ 
craft.  These  early  developments  were  the  result  of  applying  matrix 
methods,  which  were  successful  with  discrete  structures,  to  continuous 
ones.  The  first  description  of  the  procedure  was  presented  in  1956  by  a 
group  (  Turner,  Clough,  Kartin,  and  Topp)  at  the  Boeing  Aircraft  Com¬ 
pany.  The  extension  of  the  finite  element  method  to  non-structural 
problems  began  in  the  early  1960's.  Finite  element  analysis  is  closely 
related  to  the  classical  variational  concepts  of  the  Rayleigh-Ritz 
method  or  the  weighted  residual  methods  modeled  after  the  method  of 
Galerkin,  thus  it  has  been  established  as  an  important  branch  of  approx¬ 
imation  theory  (  Ref  2  and  3) • 

The  first  paper  to  propose  the  extension  of  the  finite  element 
method  to  continuum  problems  involving  the  Laplace  equation  was  by 
Zienkiewiez  and  Cheung  in  1965  (  Ref  4).  In  1968,  Martin  (  Ref  5)  using 
linear  triangular  elements  and  a  variational  principle  formulation 
solved  for  the  stream  function  for  incompressible  flow  over  a  circular 
cylinder  between  two  parallel  walls.  Norrie  and  de  Vries  (  Ref  6-8)  in 
I969  developed  finite  element  techniques  for  solving  incompressible 
flow  problems  over  single  and  multiple  airfoils.  They  used  a  variational 
principle  to  produce  the  finite  element  equations  and  formulated  the 
problem  in  terms  of  both  velocity  potential  and  stream  functions  with 
linear  triangular  elements.  The  velocity  potential  solution  involved  the 
linear  superposition  of  a  thickness  problem  and  a  lifting  problem.  The 
combination  of  the  two  problems  and  application  of  the  Kutta  condition 
at  the  sharp  trailing  edge  resulted  in  the  specification  of  the  circu¬ 
lation.  This  approach  was  also  used  by  Carey  (  Ref  9)  who  extended  it  to 
compressible  flow.. 


Shen  (  Ref  10)  formulated  the  problem  of  incompressible  flow  over  a 
lifting  airfoil  in  terms  of  the  stream  function  using  a  variational 
principle  to  develop  the  finite  element  equations.  In  this  case,  Shen 
treated  the  infinite  domain  as  an  inner  and  outer  patch.  The  inner  patch 
contained  the  airfoil  and  a  portion  of  the  flowfield  at  an  arbitrary  but 
sufficiently  distant  boundary.  Only  this  inner  patch  was  solved  through 
finite  element  procedures,  with  linear  triangular  elements.  In  the  outer 
patch  an  analytic  solution  with  unknown  coefficients  was  used.  The  two 
solutions  were  matched  at  the  common  boundary.  The  arbitrary  airfoil  was 
transformed  through  a  Joukowski  transformation  to  a  near  circle.  This 
procedure  was  also  used  by  Habashi  (Ref  ll)  who  used  a  Laurent  series 
for  the  solution  in  the  outer  patch.  He  notes  that  this  technique  pro¬ 
duces  the  value  of  circulation  without  integrating  the  pressure  distri¬ 
bution,  and  serves  to  magnify  regions  of  high  gradients,  at  the  leading 
and  trailing  edges.  Recently,  Baskharone  and  Hamed  (  Ref  12)  developed  a 
technique  for  treating  the  circulation  around  the  airfoil  as  an  additional 
variable,  the  value  of  which  is  directly  calculated  in  the  solution 
process,  rather  than  an  externally  imposed  condition.  This  is  done  by 
treating' the  circulation  as  a  nodeless  variable.  A  potential  function 
formulation  was  used  with  linear  triangular  elements  and  parabolic 
quadrilateral  isoparametric  elements.  Results  are  given  for  both  a 
single  airfoil  and  a  cascade  airfoil. 


Objective 


The  purpose  of  this  work  was  to  compare  the  ability  of  several 


different  finite  element  formulations  to  predict  the  surface  pressure 
distributions  for  two-dimensional  potential  flow  over  a  thin  airfoil  in 
an  infinite  uniform  flowfield  for  incompressible  subsonic  flow.  The 


results  were  compared  with  known  exact  solutions.  The  assumptions 
traditionally  used  in  classical  thin-airfoil  and  small  disturbance 
theories  were  used  wherever  necessary  to  simplify  the  problem.  Three 
types  of  rectangular  finite  elements  were  usedt  linear;  mixed,  a  tran 
sition  element  with  three  linear  sides  and  one  quadratic  side  aligned 
along  the  airfoil  surface;  and  quadratic,  la  grange,  elements. 
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II  Thin  Airfoil  Problem  Formulation 


The  problem  being  considered  is  one  of  steady,  two-dimensional, 
incompressible,  inviscid,  irrotationai  flow  about  a  thin  airfoil  in  a 
uniform  stream  of  infinite  domain.  Let  the  freestream  be  taken  to  be 
directed  in  the  positive  x-direction,  with  the  coordinate  system  at¬ 
tached  at  the  midpoint  of  the  airfoil.  The  symbol  n  denotes  the  infin¬ 
ite  flowfield  domain,  with  points  (x,y),  as  shown  in  Fig  1.  The  boundary 
of  J~L  (denoted  )  is  composed  of  all  points  on  the  airfoil  surface 
and  the  boundary  at  infinity.  Because  of  circulation,  a  branch  cut  is 
placed  in  .fl.  . 

For  an  irrotationai  flowfield,  a  velocity  potential  exists,  such 
that  the  governing  differential  equation  for  incompressible  flow  is  the 
Laplace  equation 

V  -  o  for  (x,y)  in  jQ_  (l) 

where  ^  is  the  velocity  potential  function.  The  boundary  conditions  are 

V§  •  VF  r  O  5  -  O  (2) 

where  R*^V-o  describes  the  airfoil  profile,  and  the  infinity  condition 

V<j)  — ?  (u*,  i  as  (x,y)  — ?  00 

where  UL  is  the  free-stream  velocity.  For  the  lifting  case,  the  Kutta 
condition  must  be  satisfied,  which  requires  the  circulation  (  P  )  be 
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such  that  the  velocity  is  finite  at  the  trailing  edge. 

Method  of  Small  Perturbations 

The  total  velocity  potential  can  be  expressed  as  the  sum  of  the 
velocity  potential  due  to  the  free  stream  and  a  perturbation  potential  0 
(Ref  13), 


(j)  -  UL( *  +  0") 


where  IL  is  the  free  stream  velocity  and  0  is  the  perturbation  poten¬ 
tial.  The  perturbation  potential  has  velocity  components  (X  and  ir 
in  the  x-  and  y-  direction  respectively  defined  by 


-  30 


LT  - 


Taking  partial  derivatives 


=  U.(  i  +  U.0, 


§_  -  U.  & 


Since  =  ixx  +  &*,  -  O  then  also 

V  0  -  o  for  (x,y)  in 


The  boundary  condition  on  the  airfoil  surface  is  now 


(  IL  +  v  V  p  r  o  i  Rx,^  =  O  (8) 


t 


The  dot  product  is  then 


(  [  \  4  "\iLE  4  -  O  (  \ 

\  2x  /  <3X  4  3Q  '  °  (9) 

The  airfoil  surface  is  described  by 

~  ^60  (10) 

then  "  f^OO  -  ^  '  O  (li) 

and  3JF  _  c(^  5F  _  _  )  (12) 

<2X  c/x 

Substituting  (12)  into  (9)  results  in 


Making  the  assumption  that  the  velocity  perturbations  produced  by  the 
airfoil  in  the  flowing  stream  are  small,  since  the  airfoil  is  assumed 
very  thin,  ^  can  be  neglected.  As  in  classical  thin -airfoil  theory, 
the  tangency  (  or  surface)  boundary  condition  is  applied  on  the  x-axis 
(y=0),  rather  than  on  the  actual  surface  itself.  The  infinity  condition 
requires 


<30  — >  o 


as  x,y  — ?  (14) 


The  Kutta  condition  requires  the  velocity  to  be  finite  at  the  trailing 
edge.  This  is  accortplished  by  setting 


at  the  trailing  edge  (15) 


In  order  to  keep  (fa  unique,  a  branch  cut  is  introduced  into  the 
flowfield.  The  potential  jump  across  this  cut  is  equal  to  the  circulation, 
defined  by 


(16) 


where  V  is  the  tangent  velocity  and  ds  is  along  the  path  of  integra¬ 
tion.  Therefore  across  the  cut 


P  --  -  0' 


is  enforced. 


(17) 


The  actual  value  of  P  is  determined  by  enforcing  the  Kutta  condition. 
The  pressure  coefficient  from  small  perturbation  theory  is  (Ref  13) 


(18) 


Due  to  the  linearity  of  the  Laplace  equation,  the  solution  for  a 
particular  airfoil  shape  can  be  determined  from  the  sum  of  the  solutions 
to  the  three  subproblems;  thickness,  camber,  and  flat  plate  at  angle  of 
attack,  as  shown  in  Fig  2.  The  difference  among  the  three  problems  is 
the  treatment  of  the  boundary  conditio~i  on  the  airfoil  (Ref  14).  For  the 
thickness  problem 


(19) 
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where 


%'?(%- It) 


and  1?^  =  upper  surface 

=  lower  surface 
For  the  camber  problem 


where 


The  condition  for  a  flat  plate  at  angle  of  attack  p(  is 


(20) 


(21) 

(22) 


^4 


(23) 


The  Finite  Element  Approximation 

The  method  of  weighted  residuals  is  used  to  obtain  the  finite 
element  equations  and  an  approximate  solution  to  the  Laplace  equation 
(Ref  3) •  To  apply  the  method  is  approximated  by 

^  u 

d>  ^  0;  (24) 

f  *  U 1 

where  0 £  are  assumed  independent  functions  chosen  such  that  all  global 
boundary  conditions  are  satisfied,  ^  are  the  unknown  nodal  parameters, 
and  K)  is  the  number  of  system  nodes.  Since  0  is  an  approximate 
solution,  when  substituted  into  the  Laplace  equation,  it  does  not 
exactly  satisfy  the  equation.  The  equation  is  thus  set  equal  to  an  error, 
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GENERAL  AIRFOIL 
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THICKNESS 


CAMBER 


CL 


ANGLE  OF  ATTACK 


Figure  2.  Airfoil  Problem  Superposition 
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£  .  Thus 


-  £ 


The  method  of  weighted  residuals  attempts  to  determine  the  •  unknowns 
in  such  a  way  that  the  error  over  the  solution  domain  is  small.  This  is 
accomplished  by  forming  a  weighted  average  of  the  error  and  specifying 
that  this  weighted  average  vanishes  over  the  solution  domain.  Thus 


ol.a  -  o 


where  w i  are  the  linearly  independent  weighting  functions.  Equation  (26) 
is  integrated  by  parts  to  introduce  the  influence  of  the  natural 
boundary  conditions 

|fv$-VW;d.a-  \  V^  nWiola-n-  (27) 

jx  a  si 

The  weight  functions  w C  are  set  equal  to  the  trial  functions 
(Galerkin's  Method).  Equation  (27)  becomes 

^  V^-VAJcol-ft  -  ^  o  (28) 


Since  equation  (28)  holds  for  any  point  in  the  domain,  it  also  holds 
for  any  group  of  points  defining  an  arbitrary  element  &  contained  in 
the  whole  domain.  Because  of  this,  a  local  approximation  can  be  made 
one  element  at  a  time,  providing  0  and  satisfy  certain  continuity 
requirements . 
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The  trial  function  for  each  element  can  be  denoted  by 


(e)  os')  C«0 

s  (29) 


Substituting  into  equation  (28) 


^  Vk)tdn.  -  ^  V(/>^  r\  K)l  d  3^  = 

SL 


(30) 


where  denotes  that  this  integral  is  non-zero  only  for  those  elements 
that  border  the  boundary  of  JT1  .  In  matrix  form,  the  element  equations 
can  be  written  as 


(31) 


where  Is  "the  elemental  stiffness  matrix,  and  is  the 

element  forcing  vector. 

iv  C  ^ 

The  elemental  K£j  and  matrices  are  described  in  Appendices  A 
through  G.  The  elemental  equations  are  transformed  into  a  global  system 
of  equations  through  an  assembly  procedure.  The  global  system  can  be 
expressed  as 


Ki-4  fi  (32) 

Treatment  of  Parfleld  Boundary 

In  the  finite  element  solution  for  the  governing  differential 
equation,  the  flowfield  is  discretized  into  a  finite  number  of  elements. 
Thus  the  infinite  domain  is  replaced  by  a  finite  domain  flp  .  There 
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are  two  possible  techniques  to  treat  the  farfield  boundary  conditions 


(Ref  15). 

The  first  method  assumes  that  nf  is  very  large,  and  that  the  actual 
gradient  boundary  conditions,  9$-?  O  ,  are  enforced  along  3X1  f  .  This 
condition  is  substituted  into  the  line  integral  term  of  the  element 
equations  ( 30 ) •  Solving  the  differential  equation  with  this  condition 
and  the  solid  boundary  condition  on  the  airfoil  (both  Neumann  boundary 
conditions)  results  in  a  non-unique  solution.  To  obtain  a  unique  solution, 
the  value  of  ^  must  be  specified  at  a  minimum  of  one  node. 

The  second  method  is  to  impose  the  condition  along  the  far¬ 

field  boundary  31^  .  The  values  of  should  ensure  satisfaction  of  the 
infinity  condition  and  be  valid  everywhere  in  the  farfield.  This  method 
has  the  advantage  of  reducing  the  solution  domain,  and  has  commonly  been 
utilized  by  others  using  finite  element  methods.  The  reduction  is  accom¬ 
plished  through  partitioning  of  the  global  system  of  equations  as  follows 


K 

k 


i 

aa  i 


Kat> 

Ku 


(33) 


The  vector  £$^is  composed  of  the  P  nodal  values  of  which  are  com¬ 
puted  from  the  expression  for  .  The  vector  is  composed  of 
the  remaining  M=N-P  unknown  nodal  values  of  ,  which  are  determined 
from  equation  (33)  by  inverting  matrix  resulting  in 


--  [W  ]  (  m  Kaw]  ( ») 


The  method  used  in  this  report  was  the  first.  The  uniqueness  of  the 


flowfield  was  established  through  specification  of  the  circulation. 


Ill  Solution  Methods 

Three  different  problem  solution  methods,  involving  different 
aspects  of  treating  the  branch  cut  were  used  to  solve  the  flat-plate 
lifting  problem. 

Superposition 

The  first  method  used  was  developed  by  de  Vries  and  Norrie  (Ref  7)* 
In  this  method,  the  velocity  potential  is  defined  as  a  superposition  of 
a  thickness  problem  and  a  lifting  problem  expressed  by 

<j>  -  P  0,  +  02  (35) 

where  V  is  the  circulation,  0  is  the  perturbation  potential,  (f> , 
and  are  the  perturbation  potentials  of  the  subproblems.  The  two 
subproblems  are  solved  as  two  separate  boundary  value  problems  stated: 
a) 


--  o 


cTn 


*  O 


<*<t> i  _ 
“cTn  ' 


O 


(0.X  • 1 


in  -TL 


on  surface  of  airfoil 


(36) 


on  farfield  boundary 

(flOu  --  ° 


where  are  the  potentials  of  the  top  branch  cut  nodes,  and  (A 

are  the  potentials  of  the  bottom  branch  cut  nodes. 
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b) 


■-  o 


in  S L 


dfl  ”  ax 


on  the  farfield  boundary 


(3?) 


on  the  surface 


(<&>„  --  o  -  o 

By  differentiating  equation  (35)  with  respect  to  x,  the  x-component 
of  velocity,  C(  ,  is  obtained. 


u  --  nut  +  Uz.  (38) 

where  U  is  the  total  x-component  velocity,  U ,  is  the  x-component 
velocity  of  (j)x  field,  and  Ui  is  the  x-component  velocity  of  ^  field. 

The  Kutta  condition  is  applied  at  the  trailing  edge  by  setting  the 
total  x-component  velocity  equal  to  zero  there.  The  circulation  can  then 
be  determined  by  equating  equation  (38)  to  zero,  thus 


pa, 


(to 


+  U-L 


(re") 


=  O 


(39) 


therefore, 


P  -  -  Lrty 

»  (TlO 


(40) 


Once  is  determined,  its  value  along  with  the  computed  values 
for  and  ^  can  be  combined  in  equation  (35)  to  determine  the  total 
perturbation  potential  field.  In  this  method,  the  velocity  along 
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the  branch  cut  is  constant.  An  advantage  of  this  method  is  that  the 


circulation  is  easily  determined  as  a  result  of  combining  (fa  and  (fa^  . 
Another  advantage  is  that  the  stiffness  matrix  for  both  and  (j>  are 
symmetric  and  banded  storage  and  solution  methods  are  applicable.  A 
Gaussian  elimination  method  was  used  to  solve  the  system  equations. 
Iterative  Method  1 

In  this  method  of  solving  the  flat-plate  problem,  the  finite 
element  equations  were  modified  for  the  jump  in  across  the  branch 
cut  by  setting 


-  4  n 

0u-  -4^ 


(41) 


and  iterating  the  solution  through  the  value  of  circulation  until  the 
U  velocity  at  the  trailing  edge  becomes  zero.  This  method  treats 
the  flowfield  as  one  boundary  value  problem,  and  requires  solution  of 
the  finite  element  equations  for  each  iteration  of  P  .  The  same  solution 
routine  as  previously  noted  was  used.  The  stiffness  matrix  is  again 
symmetric  and  banded,  thus  allowing  the  reduced  storage  method.  As  in  the 
first  method,  the  U  velocity  component  along  the  branch  cut  is 
constant . 

Iterative  Solution  Method  2 

In  this  method,  the  flowfield  is  treated  as  one  boundary  value 
problem  with  the  finite  element  equations  modified  for  the  jump  in  po¬ 
tential  <fa  across  the  branch  cut  by  defining 


h  -  A.  -  n 


(42) 
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where  the  value  of  P  is  again  determined  through  an  iterative  process 
as  in  che  previous  method.  An  additional  constraint  is  added  at  the 
trailing  edge  to  make  the  average  velocity  there  equal  to  zero, 


Vu  +VL  =  O  m 

The  requirements  of  equations  (42)  and  (43)  make  the  stiffness 
matrix  non -symmetric,  and  thus  the  reduced  storage  methods  used  in  the 
previous  methods  cannot  be  used.  This  problem  was  solved  again  using  a 
Gaussian  elimination  method.  The  U  velocity  along  the  branch  cut; 
however,  is  no  longer  required  by  the  problem  to  be  zero. 


IV  The  Flat -Plate  at  Angle  of  Attack 

The  finite  element  results  are  presented  for  a  flat-plate  airfoil 
at  a  small  angle  of  attack.  The  pressure  coefficient  distribution  over 
the  airfoil  is  determined  using  the  solutions  for  the  potential  functions. 
Appendices  A  through  G  contain  the  computational  details  of  how  this  is 
done.  The  calculated  pressure  distributions  are  compared  to  the  results 
from  classical  thin -airfoil  theory.  The  value  of  circulation,  obtained 
as  a  consequence  of  applying  the  Kutta  condition,  is  also  compared  to 
the  thin -airfoil  theory  value. 

Three  types  of  elements  were  used  to  obtain  the  pressure  coefficient 
distribution!  linear,  Appendix  A;  mixed,  Appendix  C;  and  quadratic, 
Appendix  B.  The  effect  of  the  farfield  boundary  location  and  refinement 
of  the  element  discretization  in  the  flowfield  domain  was  determined  for 
each  element.  The  superposition  method  was  used  to  generate  the  results 
in  this  section;  however,  the  iterative  solution  methods  were  shown  to 
produce  the  same  results. 

Flowfield  Grid  Parameters 

An  automatic  mesh  generation  technique  was  used  to  produce  the 
nodal  coordinates.  The  parameters  of  the  grids  were  as  follows! 

YMAX  -  the  y-direction  width  of  the  rectangular  field,  with 
the  airfoil  approximated  at  y=0. 

XMAX  -  the  x-direction  width  of  the  rectangular  field. 

NDX  -  the  number  of  divisions  in  the  x-direction  in  the  farfield, 
NDX/2  in  front  of  and  behind  the  airfoil.  The  value  of  NDX 
is  equivalent  to  one  element  for  linear  and  half  an 
element  for  quadratic  elements. 
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NDXA 


the  number  of  x-direction  divisions  along  the  airfoil 
surface.  The  same  relationship  with  respect  to  elements 
as  in  the  NDX  parameter. 

NDY  -  the  number  of  y-direction  divisions  in  the  flowfield. 

The  same  relationships  to  elements  as  in  the  x-direction 
parameters . 

Linear  Elements 

For  linear  elements,  the  U(  velocity  along  the  airfoil  surface  and 
element  boundary  interface  is  constant  for  each  element,  reference 
Appendix  A.  For  this  reason,  applying  the  Kutta  condition  at  the  trailing 
edge  requires  the  velocity  and  as  a  consequence  the  pressure  coefficient 
to  be  zero  in  the  last  element  on  the  top  and  bottom  of  the  airfoil. 

This  requirement  can  produce  a  significant  error  at  the  trailing  edge, 
if  the  last  element  is  large  with  respect  to  the  airfoil. 

Three  grid  patterns  were  used:  a  uniformly  spaced  grid,  Fig  3> 
a  uniform  grid  with  a  small  (.005  chord)  element  at  the  trailing  edge, 

Fig  4;  and  a  grid  in  which  the  size  of  the  elements  increases  exponen¬ 
tially  from  the  leading  and  trailing  edges,  in  both  the  x-  and  y-direc- 
tions  over  the  airfoil  and  in  the  field,  Fig  5* 

Farfield  Boundary  Location . 

As  previously  noted,  the  location  of  the  farfield  boundary  is 
characterized  by  the  parameters  XMAX  and  YMAX.  The  solution  of  the 
problem  could  be  obtained  for  any  combination  of  these  two  parameters; 
however,  it  is  desireable  to  reduce  the  computation  time  and  space  by 
selecting  the  smallest  possible  area.  The  method  used  was  that  used  by 
Marsh  (Ref  15).  The  lower  bounds  on  these  two  parameters,  XMAX  and  YMAX, 
were  determined  from  a  series  of  solutions  for  elements  of  fixed  size. 
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This  fixed  size,  or  Aspect  Ratio  (AR),  is  the  ratio  of  °f  the 

element.  This  analysis  was  accomplished  with  the  equally  spaced  grid 
with  small  trailing  edge  element.  Although  this  small  element  is  not 
the  same  size  as  the  others,  it  actually  has  a  beneficial  result,  as 
will  be  seen  later.  A  base  solution  for  relatively  small  values  of  XMAX 
and  YMAX,  with  constant  AR  for  the  elements  was  first  established.  The 
size  of  the  field  was  then  expanded  on  the  outer  perimeter  of  the  flow 
domain.  A  series  of  solutions  was  thus  obtained  for  increasing  values 
of  XMAX  and  YMAX.  By  holding  the  element  size  constant  and  -varying  only 
the  number  of  elements,  the  effect  of  the  farfield  boundary  location  on 
the  pressure  distribution  can  be  determined.  The  base  flowfield  was 
first  expanded  in  the  x-direction  only,  keeping  the  number  of  y-direction 
elements  constant.  The  results  of  this  expansion  are  shown  in  Fig  6. 

These  results  indicate  that  there  is  no  change  in  the  circulation  for 
XMAX  greater  than  2.5  chord! engths .  The  flowfield  was  next  expanded  in 
the  y-direction,  holding  XMAX=  2.5c.  The  results  of  this  expansion  are 
shown  in  Fig  ?.  It  is  seen  from  these  results  that  a  YMAX  greater  than 
3.0c  produces  no  change  in  the  circulation.  It  can;  therefore  be  con¬ 
cluded  that  a  farfield  boundary  location  greater  than  XMAX=  2>5c  and 
YMAX=  3«0c  has  no  effect  on  the  pressure  distribution  on  the  airfoil. 

Pressure  Coefficient ■ 

The  effect  of  the  three  types  of  grids  on  the  pressure  distribution 
was  investigated  next.  The  value  of  the  circulation  being  proportional 
to  the  integrated  pressure  distribution  can  be  used  as  a  measure  of  the 
relative  accuracy  of  the  pressure  distribution,  in  an  average  sense. 

In  Fig  8,  the  three  grids  are  compared  with  respect  to  the  convergence 
of  the  circulation ■ to  the  classical  thin-airfoil  theory  value,  when 


24 


P  (Circulation) 


Plat  Plate  at  Angle  of  Attack 
Linear  Elements  w  TE  Correction 
NDXA=  8 
NDY=  10 

T.A,T.=  thin-airfoil  theory 


T 


1  ‘  "r  —  »  —  i  -  i  i  i 

2  3  4 

XMAX 


Symbol 

AR 

YMAX 

A 

1.00 

1.25 

© 

0.50 

2.5 

□ 

0.25 

5.0 

Figure  6.  Linear  Elements,  Constant  Aspect  Ratio,  Expansion  of  Domain 
in  x-direction 
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a  uniform  reduction  in  grid  size  is  made.  Prom  this  figure,  it  is  seen 
that  the  equally  spaced  grid  is  the  slowest  of  the  three  to  converge.  An 
examination  of  the  actual  pressure  distribution  in  Fig  9  shows  that  the 
zero  pressure  requirement  in  the  last  element,  to  impose  the  Kutta 
condition,  effectively  shortens  the  airfoil  by  the  width  of  that  element. 
The  decrease  in  width  of  the  airfoil  surface  elements  can  improve  this 
condition;  however,  since  the  size  of  this  last  element  is  inversely 
proportional  to  the  number  of  elements,  it  would  take  20  elements  to 
extend  to  9 5%  of  chord  and  200  elements  to  cover  99*5%  of  chord.  From 
Fig  8,  it  is  seen  that  introducing  the  small  element  at  the  trailing 
edge  improves  the  convergence  of  the  equally  spaced  element  grid.  The 
pressure  distribution  for  this  case  is  shown  in  Fig  10.  In  this  case 
99-5^  of  the  airfoil  is  covered  by  non-zero  elements,  for  any  desired 
number  of  surface  elements.  Comparing  Figs  9  and  10,  it  can  be  seen  that 
since  the  grid  with  the  trailing  edge  element  is  effectively  longer,  for 
the  same  grid  parameters,  all  values  of  the  pressure  distribution  are 
slightly  higher  in  each  interval.  The  third  grid,  exponentially  spaced, 
is  shown  in  Fig  8  to  be  the  best  of  the  three.  The  200  element  grid, 
corresponding  to  NDX  =  NDY  =  NDXA  =  10,  has  a  circulation  value  only 
2.2#  below  that  of  the  exact  solution,  while  the  first  two  grids  with 
these  same  parameters  are  19  .J/o  and  14.1#  respectively  below  the  exact. 
The  pressure  distribution  for  the  exponential  grid  case  is  shown  in 
Fig  11.  The  major  reason  for  the  improvement  in  circulation  is  the 
concentration  of  smaller  elements  near  the  leading  and  trailing  edges. 

The  small  elements  at  the  leading  edge  allow  for  a  better  approximation 
to  the  singularity,  while  at  the  trailing  edge  they  automatically  reduce 
the  size  of  the  zei'o  pressure  element.  The  large  elements  near  the 
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Figure  9.  Pressure  Distribution  For  Linear  Equally  Spiced  Elements 
and  Uniform  Reduction  of  Element  Size 
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Figure  10,  Pressure  Distribution  For  Linear  Equally  Spaced  Elements 
With  Small  Trailing  Edge  Element  and  Uniform  Reduction 
of  Element  Size 
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Figure  11.  Pressure  Distribution  For  Linear  Exponentially  Spaced  Elements 
and  Uniform  Reduction  of  Element  Size 
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midchord  do  not  cause  a  significant  error,  as  this  portion  of  the  exact 
curve  is  the  most  linear.  In  all  three  cases  the  pressure  distribution 
step  function  approximately  intersects  the  exact  solution  at  some  point 
in  the  interval.  Thus,  even  coarse  grids  are  not  unrealistic  in  approx¬ 
imating  the  pressure  distribution.  In  Pigs  9-11 t  the  effect  of  uniformly 
reducing  the  element  size  on  the  pressure  distribution  is  shown.  The 
fastest  convergence  occurs  at  the  midchord  and  the  slowest  occurs  at 
the  leading  edge  singularity.  The  solution  trend  in  all  cases  is  correct, 
in  that  they  approach  the  exact  solution  as  the  element  size  is  made 
smaller.  Refinement  of  the  elements  through  the  farfield  parameters, 

NDX  and  NDY,  are  shown  in  Figs  12-14  for  the  three  grids.  This  technique 
produces  a  fast  convergence  from  midchord  to  the  trailing  edge,  and 
slower  convergence  toward  the  leading  edge.  Refinement  of  the  elements 
through  the  airfoil  element  width  parameter  NDXA  only  is  shown  in  Figs 
15-17  for  the  three  grids.  This  refinement  exhibits  slower  convergence 
at  the  leading  edge  than  the  previous  method;  however,  this  method  shows 
better  convergence  to  the  shape  of  the  exact  curve. 

These  refinement  techniques,  suggest  the  best  method  of  refinement 
would  be  to  first  refine  the  farfield  parameters  until  the  leading  edge 
element  pressure  varies  little  with  further  refinement,  and  then  refine 
through  the  surface  element  parameter  NDXA  to  obtain  the  best  approximation 
to  the  exact  curve. 

Mixed  Elements 

In  order  to  eliminate  the  step  function  characteristic  of  the 
pressure  distribution  for  the  previous  linear  element  grids,  a  mixed 
element  was  introduced  along  the  airfoil  surface,  Fig  18,  with  the 
quadratic  side  aligned  along  the  airfoil  boundary.  Refer  to  Appendix  C 
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Figure  12.  Pressure  Distribution  For  Linear  Equally  Spaced  Elements 
and  Uniform  Reduction  of  Farfield  Parameters  Only 


AC  (y=0) 


Flat  Plate  at  Angle  of  Attack 
Linear  Elements  w  TE  Correction 


Symbol 

Nodes 

Elements 

NDX 

NDXA 

NDY 

A 

81 

52 

8 

4 

o 

175 

136 

8 

8 

□ 

234 

190 

8 

10 

Figure  13.  Pressure  Distribution  For  Linear  Equally  Spaced  Elements 
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Figure  14.  Pressure  Distribution  For  Linear  Exponentially  Spaced  Elements 
and  Uniform  Reduction  of  Farfield  Parameters  Only 
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Figure  16.  Pressure  Distribution  For  Linear  Equally  Spaced  Elements 
With  Small  Trailing  Edge  Element  and  Uniform  Reduction 
of  Airfoil  Element  Width  Parameter 
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Figure  17.  Pressure  Distribution  For  Linear  Exponentially  Spaced  Elements 
and  Uniform  Reduction  of  Airfoil  Width  Parameter 
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for  formulation  of  equations  with  this  element.  With  this  modification 
the  velocity  and  as  a  consequence  the  pressure  along  the  airfoil  surface 
is  a  linear  function  in  each  element,  rather  than  a  constant.  This 
allows  the  application  of  the  Kutta  condition  at  the  last  node  on  the 
airfoil  rather  than  the  entire  last  element  as  done  previously.  Thus  the 
mixed  element  prevents  the  shortened  airfoil  effect  that  occurs  when 
linear  elements  are  used. 

Farfield  Boundary  Location. 

The  procedures  used  for  determining  the  farfield  boundary  location 
in  the  linear  element  section  were  also  used  for  the  mixed  elements.  The 
results  of  expanding  the  field  in  the  x-direction  are  shown  in  Fig  19, 
while  the  subsequent  expansion  in  the  y-direction  is  shown  in  Fig  20.  It 
can;  therefore,  be  concluded  that  a  farfield  boundary  location  greater 
than  XMAX=  2.5c  and  YMAX=  3. 5c  has  no  effect  on  the  pressure  distribution 
on  the  airfoil. 

Pressure  Coefficient. 

The  effect  of  element  size  on  the  pressure  coefficient  distribution 
and  circulation  were  accomplished  like  that  for  linear  elements.  The 
comparison  of  the  convergence  to  the  thin-airfoil  theory  value  of  circu¬ 
lation  with  those  of  the  two  equally  spaced  grids  in  the  previous  section 
are  shown  in  Fig  21.  From  this  figure  it  is  seen  that  the  mixed  elements 
on  the  airfoil  effectively  produce  the  same  circulation  as  the  grid  with 
the  small  trailing  edge  correction  element.  Thus  the  mixed  elements 
automatically  correct  for  the  trailing  edge  error  that  was  seen  to  be 
a  problem  for  the  equally  spaced  grid  in  the  previous  section.  The 
pressure  distribution  for  the  uniform  reduction  in  element  size  is 
shown  in  Fig  22.  When  this  distribution  is  compared  with  that  for  the 
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Figure  19.  Mixed  Elements,  Constant  Aspect  Ratio,  Expansion  of  Domain 
in  x-direction 
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Figure  20.  Mixed  Elements,  Constant  Aspect  Ratio,  Expansion  of  Domain 
in  y -direction 
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Figure  21.  Convergence  of  Circulation  of  Mixed  Elements  Compared  to 
Linear  Equally  Spaced  Grids 
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Figure  22,  Pressure  Distribution  For  Mixed  Elements  and  Uniform  Reduction 
of  Element  Size 


linear  elements  with  small  trailing  edge  element  in  Fig  10,  the  reason 
for  virtually  the  same  circulation  with  similar  parameters  becomes 
apparent.  The  value  of  the  pressure  coefficient  in  each  comparable 
element  at  the  midpoint  of  the  element  is  virtually  the  same  in  each 
case.  The  advantage  to  the  mixed  elements;  however,  is  the  improved 
approximation  compared  to  the  exact  curve  that  is  obtained  with  the 
linear  distribution  within  the  element.  Another  characteristic  of  the 
mixed  element  is  the  jump  in  pressure  at  the  interface  of  each  element. 
This  jump  is  very  small  except  at  and  near  the  leading  edge.  The  cause 
of  this  jump,  is  due  to  the  large  change  in  the  slope  of  the  exact 
solution  near  the  leading  edge.  In  attempting  to  approximate  this  change, 
the  elemental  slopes  are  also  changing  by  large  amounts  in  this  region. 
This  causes  the  potentials  at  the  common  nodes  to  vary.  The  conclusions 
of  the  previous  section  concerning  the  convergence  of  the  pressure 
distribution  to  the  exact  value  from  thin-airfoil  theory  apply  in  this 
case,  with  the  mixed  elements  producing  a  better  approximation  to  the 
exact  distribution. 

Quadratic  Elements 

The  third  type  of  element  used  was  a  quadratic,  Lagrange,  element 
which  was  used  throughout  the  flowfield,  Fig  23.  The  elemental  equation 
formulation  is  demonstrated  in  Appendix  B.  As  in  the  mixed  elements  on 
the  airfoil,  the  quadratic  elements  allow  the  Kutta  condition  to  be 
applied  at  the  last  node  on  the  airfoil.  The  Li  velocity  and  pressure 
are  also  linear  along  the  airfoil  boundary. 

Farfield  Boundary  Location. 

The  same  procedures  used  for  determining  the  farfield  boundary 
location  for  the  linear  elements  were  also  used  for  the  quadratic 
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elements.  The  results  of  expanding  the  field  in  the  x-direction  are 
shown  in  Fig  24,  while  the  results  for  the  subsequent  expansion  in  the 
y-direction  is  shown  in  Fig  25.  In  this  case  a  farfield  boundary  location 
greater  than  XMAX=  3* 5c  and  YKAX=  3* 5c  has  no  effect  on  the  pressure 
distribution  on  the  airfoil. 

Pressure  Coefficient. 

The  effect  of  element  size  on  the  pressure  distribution  is  accom¬ 
plished  as  was  done  for  the  linear  elements.  The  result  for  uniform 
reduction  of  the  element  size  parameters  is  shown  in  Fig  26.  This  figure 
shows  that  element  refinement  produces  very  fast  convergence  from  the 
quarter  chord  to  the  trailing  edge.  Near  the  leading  edge  convergence 
is  much  slower,  as  the  slope  of  the  approximation  is  changing  very  fast 
in  this  region.  In  general  it  can  be  said  that  even  a  coarse  mesh  pro¬ 
duces  a  reasonable  approximation  to  the  exact  thin -airfoil  theory  value. 
As  elements  are  refined  it  is  evident  that  the  error  in  the  pressure 
becomes  smaller.  This  fact  is  also  reflected  in  the  improvement  of 
circulation.  The  trend  of  the  solution  is  seen  to  approach  the  exact  as 
element  size  is  reduced,  which  is  the  desireable  trend. 

The  ^results  for  reduction  of  the  farfield  size  parameters  only  are 
shown  in  Fig  27.  In  this  coarse  grid,  only  four  elements  on  the  airfoil, 
there  is  little  change  in  the  pressure  distribution  from  the  quarter- 
chord  to  the  trailing  edge.  Keeping  the  airfoil  elements  constant  and 
reducing  the  farfield  elements  produces  a  significant  jump  at  the  inter¬ 
face  of  the  first  two  elements.  This  jump  increases  with  decreasing  size 
of  the  farfield  elements.  This  can  be  attributed  again  to  the  slope  in 
the  first  two  elements,  as  occurred  in  the  mixed  elements. 

The  results  for  the  reduction  of  the  element  width  over  the  airfoil 
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Figure  2?.  Pressure  Distribution  For  Quadratic  Elements  and  Uniform 
Reduction  of  Farfield  Parameters  Only 


only  is  shown  in  Fig  28.  The  results  show  that  refinement  in  this 
manner  does  not  significantly  improve  the  leading  edge  approximation, 
but  does  on  the  average  improve  the  approximation  at  each  point  to  the 
exact. 

From  these  observations,  the  same  conclusions  concerning  discretiza 
tion  that  were  made  for  linear  elements  apply  to  the  quadratic  elements. 
A  comparison  of  the  pressure  distributions  obtained  with  linear  elements 
mixed  elements,  and  quadratic  elements  with  the  element  size  parameters 
the  same  is  shown  in  Fig  29.  It  is  evident  from  this  figure  that  the 
improvement  in  circulation  is  due  to  a  closer  approximation  at  the 
leading  edge  for  the  quadratic  element.  It  is  also  seen  that  when  using 
the  midpoint  of  the  linear  elements  for  comparison,  there  is  little 
difference  in  the  three  elements  from  the  quarterchord  to  the  trailing 


Flat  Plate  at  Angle  of  Attack 
Quadratic  Elements 
XMAX=YMAX=  5c 


Figure  28.  Pressure  Distribution  For  Quadratic  Elements  And  Reduction 
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V  Conclusions  and  Recommendations 

The  Galerkin  Finite  Element  Method  was  shown  to  "be  an  effective 
means  of  approximating  the  solution  to  the  steady,  two-dimensional, 
incompressible,  potential  flow  over  a  thin  lifting  airfoil.  Three  types 
of  elements  were  used  over  the  airfoil  to  investigate  the  farfield 
boundary  location  and  the  pressure  distribution  over  the  airfoil. 

The  Superposition  Method  of  solution  was  seen  to  be  the  best  of 
the  three  solution  methods  used  from  the  standpoint  of  computer  storage 
and  computation  time.  It  was  found  through  comparison  of  the  results  of 
this  method  with  those  of  Iterative  Method  2,  that  the  constraint  of  the 
x-velocity  being  zero  along  the  branch  cut  has  negligible  influence  on 
the  pressure  distribution  and  circulation. 

The  infinite  domain  of  the  problem  was  approximated  by  a  finite 
domain.  It  was  desireable  to  make  this  domain  as  small  as  possible 
without  affecting  the  results.  This  was  accomplished  by  expansion  of  a 
base  solution  of  constant  aspect  ratio  elements  in  first  the  x-direction 
and  then  the  y-direction.  The  result  of  this  procedure  showed  that  the 
infinite  domain  could  be  approximated  by  a  relatively  small  finite 
domain.  The  domain  of  the  linear  elements  could  be  approximated  by 
XMAX=2.5c  and  YMAX=3.0c.  The  domain  of  the  mixed  elements  could  be 
approximated  by  XMAX=2,5c  and  YMAX=3.0c,  while  the  quadratic  element 
domain  was  XMAX=3.5c  and  YMAX=3.5c.  Although  the  quadratic  element 
donain  is  slightly  larger,  the  number  of  elements  used  is  much  smaller, 
14  in  both  the  x-  and  y-directions  versus  20  and  24  in  the  x-  and  y- 
directions  respectively  for  the  linear  elements. 

The  use  of  linear  elements  throughout  the  flowfield  and  on  the 
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surface  of  the  airfoil  was  investigated  for  three  grid  patterns.  It  was 
seen  that  using  a  uniform  grid  over  the  domain  was  not  an  efficient 
method,  due  to  the  singularity  at  the  leading  edge  and  the  application 
of  the  Kutta  condition  at  the  trailing  edge,  which  effectively  made 
the  airfoil  smaller.  A  correction  to  this  grid  was  made  by  adding  a 
small  element  at  the  trailing  edge.  This  modification  made  a  significant 
improvement  in  the  convergence  of  the  pressure  distribution  to  the  exact, 
by  eliminating  the  loss  of  pressure  at  the  trailing  edge  with  coarse 
grids.  A  third  grid  that  varied  the  size  of  elements,  concentrating 
small  elements  at  the  leading  and  trailing  edges,  produced  the  fastest 
convergence  to  the  exact  circulation  value.  This  was  due  to  the  improve¬ 
ment  in  the  approximation  at  the  leading  edge  singularity. 

The  use  of  mixed  (transition)  elements  on  the  airfoil  surface  with 
linear  elements  everywhere  else  was  next  investigated.  This  element  was 
an  improvement  over  the  linear  element  in  approximating  the  pressure 
distribution.  This  was  due  to  the  linear  variation  in  the  mixed  element 
as  opposed  to  the  step  function  distribution  of  the  linear  element. 

This  element  also  eliminated  the  problem  of  applying  the  Kutta  condition 
as  was  experienced  at  the  trailing  edge  with  the  linear  elements. 
Although  this  element  produces  an  improved  pressure  distribution  over 
the  linear  elements,  the  circulation  is  effectively  the  same. 

The  use  of  quadratic,  Lagrange,  elements  over  the  airfoil  surface 
and  in  the  farfield  improved  the  convergence  of  the  circulation  to  the 
exact,  when  compared  to  uniformly  spaced  linear  and  mixed  elements,  with 
an  equivalent  number  of  nodes.  Some  care  must  be  exercized  when  using  a 
coarse  grid  over  the  airfoil  due  to  the  jump  in  pressure  at  the  element 
boundaries,  near  the  leading  edge.  A  major  disadvantage  of  this  element 
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is  the  large  number  of  degrees -of -freedom  associated  Kith  it. 

From  the  above  considerations,  it  can  be  concluded  that  a  variable 
geometry  grid  with  small  elements  concentrated  at  the  leading  and 
trailing  edges  provides  the  best  results  for  the  approximation  of  the 
circulation.  This  is  further  evidenced  by  comparing  the  quadratic 
element  results  with  the  linear  elements  at  an  equal  number  of  degrees- 
of -freedom.  For  165  nodes  the  circulation  value  for  the  quadratic 
elements  is  12. k%  below  the  thin -airfoil  theory  value,  while  the  linear 
exponentially  spaced  grid  is  only  h.2%  below  thin -airfoil  theory. 

The  pressure  distribution  is  best  approximated  by  a  quadratic 
interpolation,  the  mixed  or  quadratic  elements,  on  the  airfoil  surface. 
The  combination  of  the  mixed  elements  on  the  airfoil  surface  and  linear 
elements  in  the  farfield  provides  a  reasonable  approximation  to  the 
exact  thin-airfoil  theory,  without  the  disadvantage  of  the  extra 
degrees-of -freedom  required  by  the  quadratic  elements. 

Recommendations  For  Future  Work 

This  investigation  should  serve  as  a  basis  for  further  work  in  the 
application  of  the  Finite  Element  Method  to  the  thin-airfoil  problem. 
Extension  of  the  method  to  the  cambered  and  thickness  portions  of  the 
airfoil  problem  are  necessary.  An  analysis  of  this  problem  for  the 
convergence  of  the  circulation  and  pressure  distribution  as  was  accom¬ 
plished  for  the  flat-plate  at  angle  of  attack  should  be  done.  The  treat' 
ment  of  the  branch  cut  when  the  method  is  applied  to  the  thickness 
problem  needs  to  be  addressed.  In  this  case  the  circulation  is  zero  and 
the  x-dlrection  velocity  in  the  branch  cut  cannot  be  zero. 

The  treatment  of  the  farfield  boundary  requires  further  investiga¬ 
tion.  In  this  report  no  potential  values  were  defined  on  the  farfield 
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boundary.  A  method  for  defining  the  potentials  on  the  farfield  boundary 
would  reduce  the  system  of  equations  to  be  solved  through  partitioning. 

Further  investigation  of  grid  refinement  needs  to  be  made.  The 
concentration  of  small  elements  at  the  leading  and  trailing  edges  should 
be  extended  to  include  the  mixed  and  quadratic  elements.  The  effect  of 
placing  a  mid-side  node,  from  either  the  mixed  or  quadratic  element,  at 
the  leading  and  trailing  edges  should  be  reviewed. 

Other  types  of  elements,  such  as  Hermitian,  to  guarantee  continuity 
of  the  pressure  at  element  boundaries,  infinite  elements,  that  do  not 
require  the  finite  boundary,  or  special  elements  to  improve  the  approxi¬ 
mation  at  the  leading  edge  singularity  could  also  be  applied  to  this 
problem. 

The  orientation  of  the  branch  cut  in  this  report  was  restricted  to 
extension  horizontally  from  the  trailing  edge  to  the  farfield  boundary 
along  y=0.  The  effect  of  changing  this  orientation,  such  as  approximating 
it  with  a  function  whose  slope  at  the  trailing  edge  were  equal  to  the 
angle  of  attack  and  that  became  horizontal  at  some  distance  behind  the 
trailing  edge,  should  be  compared  with  the  results  obtained  here. 
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Finite  Element  Equations  For  Flow  Over 
an  Airfoil  For  a  Bilinear,  Rectangular 
Element 


Interpolation  Functions 

The  interpolation  functions  for  the  bilinear  rectangular  element 
shown  in  Fig  30  are  given  by 


(A-i) 

i(?-0(r+0 

(A-2) 
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Coordinates  (?,r)  are  the  local  nodal  coordinates . 

Elemental  Equations 

The  Finite  element  equations  obtained  from  the  governing  differential 
equation  and  written  in  elemental  form  are  expressed  by  equation  (31)  as 
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Global  coordinates  transformed  to  the  local  coordinates 


Figure  30*  Bilinear  Rectangular  Element 
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by  the  transformation  equations! 
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From  these  expressions,  the  elemental  equations  can  be  written  in  the 
form 
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The  matrix  A(j  is  dependent  only  on  the  size  of  the  element  and 
not  directly  on  any  particular  airfoil  contour.  Substituting  the  inter¬ 
polation  functions,  A-l  through  A -4,  into  A -5  results  in 
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The  vector  is  directly  dependent  on  the  airfoil  surface 
contour.  This  quantity  is  evaluated  only  for  elements  that  share  a 
common  boundary  with  the  surface.  For  all  other  elements,  this  vector 
is  zero. 

Velocity  Distribution 

The  velocity  in  element  £  is  calculated  from  the  assumed  solution 
for  the  potential  given  by  eq  24 


Prom  the  definition  of  the  potential  function 
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The  x-dlrection  velocity  on  the  top  surface,  .  u  becomes 
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Likewise  on  the  lower  surface,  becomes 
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Pressure  Distribution 

For  steady,  small-disturbance  theory  the  coefficient  of  pressure  is 
given  by  (Ref  13) 


Cp  -  -  2  U.  (A-17) 

From  thin-airfoil  theory,  the  pressure  coefficient  is  evaluated 
along  i^o+'  •  Thus  for  element  &.  which  borders  the  airfoil  surface, 

the  elemental  pressure  coefficient  becomes 

CF  ( ’Sb1  o*  V  ^ f *±0  (*-«> 

For  an  upper  surface  element,  -  /  ,  Cp  .becomes 

e.  , 

CPo  =  *“£,  (^4  ~  (*-*9) 

For  a  lower  surface  element,  ^=+|  ,  Cp  becomes 

CPc  =  "  X  (  0|  “  (A-20) 

Within  the  element  <£.  ,  the  pressure  is  a  constant  value,  which 
results  in  "jumps”  in  pressure  between  elements  along  the  airfoil 
surface.  The  pressure  distribution  along  is;  therefore,  a  step 

function  for  the  bilinear  rectangular  element. 

Boundary  Influence 

Those  elements  along  the  airfoil  surface  will  produce  non-zero 
forcing  terms  ..  The  forcing  term  is  determined  by  evaluating  eq  A-7 
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for  each  element 
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For  a  flat  plate  at  angle  of  attack  o( 


integral  to  the  local  coordinate  system  gives 
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For  elements  along  the  lower  surface  the  forcing  terms  are  determined 
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For  elements  along  the  upper  surface  the  forcing  terms  are  determined 
from 
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2ndix  B 


Finite  Element  Equations  For  Flow  Over  an 


Airfoil  For  a  Biquadratic.  Lagrange. 
Rectangular  Element 


Interpolation  Functions 

The  interpolation  functions  for  the  biquadratic,  Lagrange, 
rectangular  element  shown  in  Fig  %  are  given  by  (Ref  3) 
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Coordinates  CS.f)  are  the  local  nodal  coordinates. 
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Elemental  Equations 

The  elemental  equations  are  the  same  as  those  derived  for  the 
bilinear  rectangular  element,  eq  A -5  through  A-7, 
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The  first  matrix  is  only  dependent  on  the  size  of  the  element 
and  not  directly  on  any  particular  airfoil  contour.  Substituting  the 
interpolation  functions,  B-l  through  B-9,  into  A -5  results  in  the  matrix 
shown  in  Fig  32. 

The  vector  ^  is  directly  dependent  on  the  airfoil  surface 
contour.  This  quantity  is  evaluated  only  for  elements  that  share  a 
common  boundary  with  the  surface.  For  all  other  elements  this  quantity 
is  zero. 

Velocity  Distribution 

The  velocities  in  element  C.  are  calculated  as  in  Appendix  A, 
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Figure  32.  Elemental  Stiffness  Matrix  for  the  Quadratic,  Lagrange,  Element 


The  x-direction  velocity  on  the  top  airfoil  surface,  becomes 

U»  *  Oj  +  (  i  +-()<f>Z  -  2  %  ^5"]  (B-ll) 

Likewise  on  the  lower  surface,  ,  4  becomes 

i  C(S  +  4.N)  03  4  (  S‘TS)04  -  (B-12) 

Pressure  Distribution 

For  steady  small  disturbance  theory  the  coefficient  of  pressure 
is  given  by  (Ref  13) 


From  thin-airfoil  theory,  the  pressure  coefficient  is  evaluated 
along  o^-=  0*  ,  Thus  for  element  6  which  borders  the  airfoil  surface 

the  elemental  pressure  coefficient  becomes 

c® (x,q=  o1) -  k-  %  ± 0  (B-i-*) 

<J 

For  an  upper  surface  element,  ^=--1  ,  CP  becomes 

For  a  lower  surface  element,  .  cp  becomes 

ca  1  -  £[(  -Z^J  (b-16) 

Within  the  element  6.  the  pressure  varies  linearly.  Since  this 
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.  u 

element  is  a  C  element,  there  is  a  jump  in  pressure  from  element  to 
element . 

Boundary  Influence 

Elements  along  the  airfoil  surface  will  produce  non-zero  forcing 
terms  ^  forcing  terms  are  determined  as  in  Appendix  A,  and 

expressed  in  A-7  as 


Si  -  1  4^  0;1  d* 


(B-17) 


For  a  flat  plate  at  angle  of  attack  o(  ,  Converting  the 


integral  to  the  local  system  gives 


*  C  -ola>  K)t  I  ol  % 

An  1  J 


(B-18) 


For  elements  along  the  lower  surface  the  forcing  terms  are  determined 


T  « 


(B-19) 


For  elements  along  the  upper  surface  the  forcing  terms  are  determined 


-  i 

Sc  =-•<*  ^(S^-hJdV 


(B-20) 


Appendix  G 

Finite  Element  Equations  For  Flow  Over  an 
Airfoil  For  a  Mixed 
Rectangular  Element 


Interpolation  Functions 

The  interpolation  functions  for  the  mixed  rectangular  element 
shown  in  Fig  33  are  given  by 

N,  -  (C-D 

Wj-  ^■(i  +  V)0-10  (0-2) 

5  ( i+V Y  i+r")  (c-3) 

tJ4-- (°-*> 

Kb;  --  (o-5) 

Coordinates  («.*)  are  the  local  nodal  coordinates. 

Elemental  Equations 

The  elemental  equations  are  the  same  as  those  derived  for  the 
bilinear  rectangular  element,  eq  A -5  through  A -7 

ACj  <tj  -  (A-5) 

AtJ  •  §(£*$*£* tH&Tbyw  (A-6) 

5*  *  «->  (  «K  (A-7) 

-1(40  Y*±l 
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The  matrix  is  only  dependent  on  the  size  of  the  element  and 
not  directly  to  any  particular  airfoil  contour.  Substituting  the 
interpolation  functions,  C-l  through  C-5,  into  A -5  results  in 
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(c-6) 


The  vector  ^  is  directly  dependent  on  the  airfoil  surface 
contour.  This  quantity  is  evaluated  only  for  elements  that  share  a 
common  boundary  with  the  surface.  For  all  other  elements,  this  vector 
is  zero. 

Velocity  Distribution 

The  velocity  in  the  x-direction  on  the  airfoil  is  calculated  as  in 
Appendix  B  for  the  quadratic  element.  For  the  upper  surface,  C(  becomes 

Uu  +  C (c'7) 
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Pressure  Coefficient 


The  pressure  coefficient  on  the  airfoil  is  calculated  as  in 
Appendix  B  for  the  quadratic  element.  For  the  upper  surface, 
becomes 

cPe*  (c.6) 


Boundary  Influence 

Those  elements  fl.  along  the  airfoil  surface  will  produce  non-zero 
forcing  terms  .  Since  the  quadratic  side  of  the  element  is  positioned 
along  the  airfoil  surface,  the  boundary  influence  is  the  same  as  that 
derived  in  Appendix  B  for  the  quadratic  element. 
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